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Abstract 

An analytical treatment for the sedimentation rate of disordered suspensions 
is presented in the context of a resistance problem. From the calculation it is 
confirmed that the lubrication effect is important in contrast to the sugges- 
tion by Brady and Durlofsky (Phys.Fluids 31, 717 (1988)). The calculated 
sedimentation rate agrees well with the experimental results in all range of 

the volume fraction. 
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The sedimentation of disordered suspensions is important in both technology and lab- 
oratory [|IJ. The role of the sedimentation is relevant to the current topics of statistical 
mechanics such as fluidized beds of gas-solid or liquid-solid mixtures [||-|5| , and the density 
waves in the granular flows in vertical tubes ||. We believe the subject is a fundamental one 
in fluid mechanics [[/J. The rate of sedimentation for disordered suspensions under gravity 
has yet to be determined theoretically except for a problem for dilute spheres with hard core 
interactions at a small Reynolds number |L]||. 

Our present understanding of theoretical studies of sedimentation of monodisperse ran- 
dom suspensions can be summarized as follows. Batchelor || has calculated the sedimenta- 
tion rate in the dilute limit of hard core particles with the radius a based on the following 
assumptions: (i) The rate can be obtained from the combination of the mobility matrix of 
two particles and the two-body correlation function g e q{f) where r is the relative distance 
of particles, and (ii) the correlation function is assumed to be g eq (r) = 6{r — 2a) , where 
8(x) is the step function 6(x) = 1 for > and 8(x) = otherwise. His result at the volume 
fraction d> can be written as U(<fi)/U = 1 - 6.550 + O(0 2 ) for -> 0, where U(cj)) the 
sedimentation velocity at and Uq is the equilibrium sedimentation velocity of one particle. 
The result of Batchelor consists of two parts: one is 1 — 50 from the Rotne-Prager tensor 
which represents the effects of long-range hydrodynamic interaction, and another —1.550 
from the lubrication, the hydrodynamic repulsive, force. Extensions of this dilute theory to 
concentrated suspensions require the account of many body hydrodynamic interactions. A 
generalization ||, based on the method of O'Brien JK| predicts negative sedimentation rate 



for > 0.27. Brady and Durlofsky ]TT| have also obtained a negative sedimentation rate for 



> 0.23 when they adopt well accepted correlation function g e q{f) for concentrated suspen- 
sions. As a result, they claim that the Rotne-Prager approximation actually captures the 
correct features of sedimentation and ignored all of the contributions from the lubrication 
force. We feel, however, the statement by Brady and Durlofsky Jll] unacceptable, because 
there is no reason to ignore lubrication effects in the dilute limit ||. On the other hand, 
Beenakker and Mazur [|T|,|l3|] also calculated the sedimentation rate based on an effective 



medium approximation and multipole expansions. Although they did not present an ex- 



plicit expression of the sedimentation rate, Ladd |14j] indicated that their result is better 



than the result by Brady and Durlofsky [II] for concentrated suspensions. In this Rapid 
Communication, we wish to demonstrate the relevance of the lubrication force and improve 
the theory by Brady and Durlofsky [|TT|] . We also clarify the relationship between our theory 



and that by Beenakker and Mazur [|TT],|T2l . 

The problem of sedimentation of N particles with the radius a at low Reynolds numbers 
is equivalent to obtaining the resistance matrix R or the mobility matrix M in 

U = ^^M-F, M = R _1 , (1) 

where U and F denote the sets of the velocity field of N particles and the force exerted on iV 
particles, respectively, and /i is the shear viscosity. These mobility and resistance problems 
are not easy to solve even numerically. One of the most successful numerical methods, the 
Stokesian dynamics, has been developed by Brady and his coworkers fl5| , |l6]| . The extension 
by Ladd |TJ| also follows a similar algorithm to the Stokesian dynamics. They decouple the 
resistance matrix into the far-field part (M°°) _1 and the lubrication part R lub as 

R = (M°°) _1 + R lub , (2) 

where R lub is calculated by the pairwise additive expression of the two-body lubrication 
matrix R^ = R 2 b — (M^) _1 . The resistance matrix is calculated as a function of the particle 
configuration at each numerical step. Then the force exerted on spheres and consequently 
the equation of motion are obtained. The success in the Stokesian dynamics suggests that 
the problem for sedimentations should be considered based on a resistance picture. In fact, 
some unphysical results of simulations based on a mobility picture supports this statement. 
We may understand the relevance of a resistance picture as follows. Since the contribution 
of the lubrication is proportional to the number of particles, as will be shown, the direct 
addition of the lubrication for the mobility cannot avoid a negative sedimentation rate. In 
other words, the linear contribution of the lubrication to the drag is reasonable, while the 



linear addition of the lubrication to the mobility cannot produce any nonlinear complicated 
motion of particles in experiments. 

Thus, we are not surprized by the failure of direct generalizations of Batchelor's theory 
which is described as a mobility problem. We must calculate the sedimentation rate in the 
context of a resistance problem. The problem is, thus, reduced to obtaining < (M°°) _1 > 
_l_ < p/u£> > where the bracket is the average over the particle configurations. Note that 
< (M°°) _1 > and < R lub > are the scalar quantities. The far-field part can be calculated 
from < (M 00 )- 1 >~< M°° >-*= M(k = 0) _1 , where M(k) is defined by 

M{k) = 1 + n j e jkri2 (^(r 12 ) - l)k • G(r 12 ) • kdr 12 . (3) 

Here k = k/fc, the relative position r 12 of the particles 1 and 2, n is the number density of 
partcles. The explicit representation of the Fourier component of the tensor G = {G^} is 
given by HO 

G y (k) = tea*^{5« - ^) (4) 

with the spherical Bessel function j (ka). For later discussion we drop the suffix of r 12 and 
assume the isotropy of systems as g eq {r = |r|). 

The correlation function g eq {r) can be approximated |T^JT7|] by the equilibrium distribu- 
tion function for hard sphere systems based on the Percus-Yevick approximation fl8 l. The 
Fourier transform of g e q{r) — 1, h(k) is represented by [|H 



_ 47ra 3 ~c(ka) 
h{k) ~ l + 30c(te)' (5) 

where c(x) is the direct correlation function which also depends on <fi. The correlation 

function in (||) reduces to g eq (r) = 6{r — 2a) in the dilute limit. From (|5|) we can evaluate 

< M°° >= | / °°rfx(^) 2 (l + 30c(x)) _1 numerically. Brady and Durlofsky (ll] evaluated 



this ]20| by using the Laplace transform of the Percus-Yevick distribution function 0] and 



the method of O'Brien |T(J as 



4 



which is a correct evaluation of the contribution from the far-field part. 

Now, we evaluate the contribution from < R lub >. For simplicity of the argument, we 
neglect contributions from higher order moments such as torque and shear. Since < R lub > 
is evaluated from a pairwise additive approximation, < R lub > is represented by 



< R 



l Lib 



> n l^drg eq (r)k 



A 11 +A 12 -{(M^)- 1 11 + (M^)- 1 12 } 



•k. 



(7) 



The tensor A Q/ 3 is a part of R2B and its suflicies represent the particles. The tensor An + Ai 2 , 
thus, is given by 

fY u + Y u \ 



An + Ai 2 



Y n + Y 12 
V X n + X 12 J 



where the explicit representations of X« and K,- are given by Jeffrey and Onishi [22 



(8) 



as a 



series expression. On the other hand, (M^)n is the unit tensor and (M 2*5) 12 is the Rotne- 
Prager tensor which is represented by 



:MS)i2=x 00 (r)rr + 2/ 00 (r)(l-fr), 



(9) 



where x°°(r) = (3/2)(r/a)~ 1 -(r/a)~ 3 and y°°(r) = (3/4)(r/a)- 1 + (l/2)(r/a)~ 3 . The tensor 



'M 



00 \— 1 

2B) 



u + (M 



00 \— 1 

2BJ 12 



can be readily calculated as 



:M 2 °° b )- 1 11 + (M? B )-\ 2 = X°°(r)rr + Y°°(r) (I - rr) , 



(10) 



where X°°(r) = (1 + x°°(r)) 1 and Y°°(r) = (1 + y°°(r)) 1 . Thus, the average of the 
contribution from the lubrication part is described by 



poo 

< R"* >= / dz z 2 g eq (r)W(z) 

J2 



(11) 



where z = r/a and 

W(z) = X u + Xi 2 + 2F n + 2Fi 2 - 



6z 3 (-2 + 5z 2 + Az 3 



-2 + 3z 2 + 2z 3 )(2 + 3z 2 + Az 3 ) ' 
With the aid of the exact result by Jeffrey and Onishi |22| W(z) can be evaluated as 



(12) 



rzr , \ 21 1 789 1 n/ L 

•"'''"F-m? 40 '?' (13) 

for 2 ^> 1. Thus we can evaluate < R lub > by the numerical integral. For the practical 

purpose, it is convenient to have an explicit expression for < R lub >. If we assume g e q{f) = 

9{r — 2a), < R lub > is approximately represented by 

r 20 roo /21 1 789 1 \ 

< R lub >~ / dzz 2 W(z) + / dzz 2 t ~ 1.4920 (14) 

J 2 J20 \ 4 z 4 64 z 5 / 

When we compare the result ([14]) with the one obtained with the aid of the Percus-Yevick 
distribution function for g eq (r), we find that the two results have no significant difference 
(see Fig.l). This statement is applicable to the calculation for the lubrication part of the 
mobility matrix as < M'" 6 >~ —1.550. We thus confirm that the contribution from the 
lubrication is insensitive to the form of g eq {r) and is proportional to 0. Thus we should solve 
the problem in the context of a resistance problem to avoid a negative sedimentation rate. 
From (^) and flT4] ) we obtain 

^ c 1 = a-^) 3 fl5 ) 

U Q < M°° >-! + < R lub > l + 20+ 1.4920(l-0) 3 ' v ; 

As will be shown, this result is sufficiently close to experimental values. The dilute limit 
of our result U/Uq = 1 — 6.490 + O(0 2 ) is slightly different from Batchelor's result U/Uq = 
1 — 6. 550 + (0 2 ). This discrepancy comes from the relation R lub ^ (M iub ) _1 . The true dilute 
limit should be calculated under the considerations of all of higher order moments P5[ . It is 
worthwhile, however, to indicate that our theory essentially resolves the contradiction about 
contributions from the lubrication in the result by Brady and Durlofsky JCTJ]. 

Now we compare our result with that by Beenakker and Mazur fl2fl . They rewrite the 
renormalized (4) as 

M 10 (k) = l + k-G 7o (r = 0)-k + n J dr e ik ' r k ■ G 70 (r) • k {^(r) - 1} , (16) 
where G 70 (r) is given by 



dk ik-r <pS yo (ka) 
(2tt) 3 l + (j)S^(ka) 



G 70 ( r ) = G(r) - ^e^ , ^iTlL m- (17) 
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Here S l0 (x) is the structure factor and G(r) = for r = and G(r) = G(r) for 
r 7^ 0. Substituting @ into (|16|) and noting < M >= M l0 {k = 0), we obtain 
< M >= l^dx (^) 2 {(1 + cj)S J0 (x))(l + 3<f)c(x))}- 1 . The function S l0 {x) tends to 



5/2 for dilute case and small x. In the dilute limit, the result by Beenakker and Mazur (12 
is reduced to U/Uq ~ 1 — (15/2)0 + O(0 2 ), which is considerably away from Batchelor's re- 
sult ||. Even in concentrated cases, S J0 (x) still may be replaced by 5/2, although its actual 
expression is complicated. With the aid of (|6|) an approximate expression of Beenakker and 
Mazur [12| is given by 

U (1-0) 3 



U (1 + 20)(1 + 50/2)' 



(18) 



From (jig), it is easy to understand that Beenakker and Mazur |12| renormalize the Rotne- 
Prager tensor by taking into account the contribution from the structure factor. The de- 
viation from Batchelor's result in the dilute limit suggests that they miss the quantitative 
description for the short-range force, because their effective field approximation includes 
only parts of the lubrication by a collection of ladder diagrams. Their theory, however, 
may be good for dense suspensions where the requirement for their approximation may be 
satisfied. 

Let us compare theoretical results with experimental ones [p^ -p8|(Fig.2). We recognize 



that our theory improves the result by Brady and Durlofsky |TT|| and achieves good agreement 
with experiments. Therefore, we conclude that the contribution from the lubrication force 
is small but relevant. For < 0.2, it seems that our result is better than that by Beenakker 
and Mazur [12]. In high concentration regions, however, our sedimentation rate is a little 



larger than the experimental values, while the prediction by Beenakker and Mazur |T2[ works 
well. This disagreement between our theory and experiments in concentrated regions seems 
to come from the neglect of higher order moments. The high sedimentation rate without 



higher order moments for a regular configuration of particles has been reported ||15|| . To 
check this tendency for random particle configurations we have performed a simulation for 
50 particles based on the Stokesian dynamics, where we neglect the contributions from higher 



order multipole expansions. In our simulation the particle configuration is at random and 
average 100 configurations for each <fi to calculate the sedimentation rate. When we neglect 
the statistical error, the tendency of high sedimentation rates in large <ft coincides with that 
of our theory. 

In conclusion, we have confirmed that the calculation of sedimentation rate should be 
performed in the context of a resistance problem. It is not surprising that the direct gener- 
alization of Batchelor's theory based on the mobility picture gives us wrong answers. Thus, 
we should include the lubrication effects in contrast to the claim by Brady and Durlofsky 
1 Our method including the lubrication force is an adequate systematic approach to 
extend the dilute theory. We demonstrate that the lowest order contribution to the sedi- 
mentation rate of the lubrication force becomes closer to experimental values than that by 
the Rotne-Prager approximation. The discrepancy between our calculation and experiments 
at high should be improved if we include the contribution from torque and other moments. 
The consecutive improvement of our calculation of the sedimentation rate will be reported 
elsewhere. 
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FIGURES 

FIG. 1. The comparison of several theoretical and numerical predictions of the sedimentation 
rate U{4>)/Uq as functions of eft. For Eq.(ll) with PY, we use the Percus-Yevick distribution 
function for g e q{r) in (11) to evaluate < R lub >. The result of Ref. (T^j is obtained from his precise 
simulation. Ref. Jl^] is from their Fig. 2 with k = and its approximate expression is given by (|l8|). 

FIG. 2. The comparison of several theoretical results with experimental results for U((J))/Uq. 
We also plot the data of our Monte-Carlo simulation. See the text for the details. 
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